package theory.vanpelt.bifurcation;

import common.statistics.Series;

/**
 * Estimator for the probability density of BE-model.
 * 
 * The density is computed by recursive function. To save computation time, 
 * the simulation time is set to 100 steps, and trees with degree larger 
 * than 100 are ignored. Note that this class uses Comb and Fact for 
 * pre-computed combinatorial / factor values. The pre-computation goes 
 * up to 100. If the range of values is changed here, then Comb/Fact have
 * to be changed accordingly.
 *   
 * Ref: Van Pelt papers (year 200x)
 * 
 * @author risto
 *
 */
public class PDistEstimate {

	double B = 10, E = 0.6;
	public int T = 100;
	public int MAX_DEG = 100;

	public double P[][] = new double[MAX_DEG][T+1];
	public double EX[] = new double[T+1];
	public double dev[] = new double[T+1];

	double p(int n)
	{
		return B/(T*Math.pow(n, E));
	}
/*
	BigInteger fact(int a)
	{
            if (a < 2) return BigInteger.ONE;

            BigInteger n = BigInteger.ONE;
            for (int i=1; i<=a; i++) {
                n = n.multiply(BigInteger.valueOf(i));
            }
            return n;
	}
*/

	public PDistEstimate(double B, double E) {

		this.B = B;
		this.E = E;

		// init; the initial tree has P(degree=1) = 1, P(degree!=1) = 0.
		for (int i = 0; i < MAX_DEG; i++) P[i][1] = 0;
		P[1][1] = 1;

		double fperror = 0;

		// loop; for each generation, compute the probabilities for each P(degree=n).
		for (int i = 2; i <= T; i++) { // generation
			
			double esum = 0;
			double devsum = 0;
			
			for (int n = 1; n < MAX_DEG; n++) { // tree degree
				
				// P(n,i) = sum_j=0:n/2 [ P(n-j,i-1) * (n-j  j) * p(n-j)^j * (1-p(n-j))^(n-2j) ]
				double sum = 0;
				//int jMax = (int) Math.ceil(n/1.0);
				//for (int j = 0; j < jMax; j++) {
                for (int j = 0; j < n; j++) { // in van Pelt paper this was n/2, but that didn't seem to work??
					double Pprev = P[n-j][i-1];
					double c = Comb.comb(n-j, j).doubleValue();
//                    double c = Comb.comb[(n-j) + j*100].doubleValue();
//                    System.out.println(c);
                    double pnj = B/(T*Math.pow(n-j, E));
//                    System.out.println(pnj);
					double pa = Math.pow(pnj,j);
					double pb = Math.pow(1-pnj, n-2*j);
//					double pa = Math.pow(p(n-j),j);
//					double pb = Math.pow(1-p(n-j), n-2*j);
//					double sumDelta = Pprev * c * pa * pb; 
	//				double foo = Math.log(Pprev) + Math.log(c) + Math.log(pa) + Math.log(pb);
//					double diff = Math.exp(foo) - sumDelta;
		//			fperror += Math.abs(diff);
//					if (diff > 0) System.out.format("%.50f - %.50f\n", sumDelta, diff);
					sum += Pprev * c * pa * pb;
//					sum += Math.exp(foo);
				}
                
                // probability of tree of degree n in i:th generation
				P[n][i] = sum;
				
				// early stop
				if (sum <= 0.0001) break;
				
//				if ((n > MAX_DEG*2/4) && (sum < 0.0001)) break;
				
				//				System.out.println("P[n:"+n+"][i:"+i+"] = " + sum);
				esum += sum * n;
			}
			
			// expected degree of tree
			EX[i] = esum;

			// compute deviation
			for (int n = 1; n < MAX_DEG; n++) {
				double pn = P[n][i];
				devsum += Math.pow(n-esum,2)*pn;
			}
			dev[i] = Math.sqrt(devsum);

			//System.out.println("E["+i+"] = " + esum + ",  +- " + dev[i]);

		}

//		System.out.println("Floating point error sum " + fperror);
		
	}

	@Override
	public String toString() {
		return String.format("PDist(%.4f,%.4f)->(%.4f,%.4f)", B, E, getMeanEst(), getDevEst());
	}
	
	public double getMeanEst()
	{
		return EX[EX.length-1];
	}

	public double getDevEst()
	{
		return dev[dev.length-1];
	}
	
	public double sample()
	{
		//P[n][i] = sum;   n= 1 .. MAXDEG-1,  i=T
		double sum = 0;
		for (int i = 1; i < MAX_DEG; i++) sum += P[i][T];
		
		double X = Math.random();
		double cdf = 0;
		for (int i = 1; i < MAX_DEG; i++) {
			cdf += P[i][T]/sum;
			if (X <= cdf) return i;
		}
		return MAX_DEG;
	}
	
	public Series sample(int n)
	{
		Series s = new Series();
		while(n-- > 0) s.add(sample());
		return s;
	}


	public static void main(String[] args) throws Exception {
		
		long start = System.currentTimeMillis();
		PDistEstimate pd = new PDistEstimate(8, 0.2);
		System.out.println("Mean " + pd.getMeanEst());
		long end = System.currentTimeMillis();
		System.out.println("TIME " + (end-start));
		
		/*
		// compute distributions for various parameter combinations, and save results to pdist.txt
		FileWriter fw = new FileWriter("pdist.txt");
		for (int iE = 0; iE < 100; iE += 2) {
			for (int iB = 15; iB < 60; iB += 1) {
				double B = iB/10.0;
				double E = iE/100.0;
				PDistEstimate pd = new PDistEstimate(B, E);
				String out = String.format("%.2f,%.2f=%.2f,%.2f\n", B, E, pd.getMeanEst(), pd.getDevEst());
				System.out.print(out);
				fw.write(out);
			}
			fw.flush();
		}
		fw.flush();
		fw.close();
		*/
	}
}
